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We investigate the linear stability of a flat interface that separates a liquid layer 
from a fully-developed turbulent gas flow. In this context, linear-stability analysis 
involves the study of the dynamics of a small-amplitude wave on the interface, and we 
develop a model that describes wave-induced perturbation turbulent stresses (PTS). 
We demonstrate the effect of the PTS on the stability properties of the system in 
two cases: for a laminar thin film, and for deep-water waves. In the first case, we 
find that the PTS have little effect on the growth rate of the waves, although they 
do affect the structure of the perturbation velocities. In the second case, the PTS 
enhance the maximum growth rate, although the overall shape of the dispersion curve 
is unchanged. Again, the PTS modify the structure of the velocity field, especially 
at longer wavelengths. Finally, we demonstrate a kind of parameter tuning that 
enables the production of the thin-film (slow) waves in a deep-water setting. 

I. INTRODUCTION 

Understanding the stability of a flat interface separating the liquid and gas phases in a 
turbulent boundary layer is a long-standing problem [T]. In this paper, we focus in detail on 
one aspect of this problem, namely modelling the perturbation turbulent stresses that arise 
when a wave on the interface interacts with the turbulence and generates Reynolds stresses. 
The approach we take has multiple facets: we derive a model velocity field describing the 
base turbulent flow in the system, and we then perform a linear-stability analysis around 
this base state using the Orr^Sommerfeld formalism. We construct a detailed model to 
describe the turbulent stresses that appear in these linear-stability equations. Motivated by 
the paper of lerley and Miles [2], we use a model that linearly couples the streamfunction 
to the turbulent kinetic energy and Reynolds stresses, and divides the gas domain into 
near-equilibrium and rapid-distortion regions. 

The purpose of our investigation is twofold. In the present paper, we develop a model 
that describes in detail the averaged effects of turbulence on the linear stability of the gas- 
liquid interface. We carry out some preliminary calculations for two archetypal cases: the 
thin liquid film sheared by a turbulent gas flow, and the generation of interfacial waves over 
deep water. This model forms the basis for a more detailed study of the thin- film case (in 
a pressure-driven scenario), which is reported elsewhere [3j. Since the present work focusses 
on turbulent gas-liquid instability in general, we place our work in context by describing the 
studies that have been performed on similar gas-liquid systems in the past. 
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The generation of waves by wind: The problem of computing the growth rate of an 
interfacial instabihty for a turbulent gas-liquid system was first considered by Jeffreys [4J, 
Phillips [i5j , and Miles [6j . Phillips studied a resonant interaction between turbulent pressure 
fluctuations and capillary- gravity waves on an interface, while Miles focussed on a particular 
kind of shear inviscid flow that produces linear instability. In the paper of Miles, in later 
works by the same author 0|8j, and in the work of Benjamin [9j, the liquid layer is neglected, 
and is replaced by a wavy wall, and the problem is reduced to computing the stress exerted 
on the interface by the turbulent gas flow. Energy transfer from the gas to the interface is 
governed by the second derivative of the mean flow. Indeed, the growth rate of the instability 
is determined by the sign of the second derivative of the mean flow at the critical height - 
the height at which the mean flow and wave speed are equal. Recent work of Lin et al [10] 
indicates that the Phillips mechanism may be important in an initial regime, leading to 
wave growth that is approximately linear in time. Their results show that this is followed 
by an exponential growth regime, which is primarily governed by the disturbances in the 
flow induced by the waves themselves. In this paper, we focus on this latter stage; we shall 
investigate the effects of the presence of short waves (possibly induced by a Phillips-type 
mechanism) on this exponential wave growth regime in a later, more detailed, study. 

Modelling Reynolds stresses: In these early works O [71 El |9] , the turbulent nature of 
the gas flow is taken into account through the prescription of a logarithmic mean proflle in 
the gas. Nevertheless, the Reynolds stress terms that enter into the stability equations are 
ignored. This problem is rectifled by Van Duin and Janssen [llj, and by Belcher and co- 
workers in a series of papers J2l[13l[ll]. The latter group studies the interfacial stability of 
a sheared air-water interface, and specialize to an air-water system for oceanographical ap- 
plications. With the exception of [T3], they treat the interface in a manner similar to Miles. 
Particular care is taken in developing an understanding of the structure of the turbulent 
shear stresses inherent in the problem through the use of scaling arguments and a truncated 
mixing-length model. This is representative of an approximation of a Reynolds- averaged en- 
semble of realizations of the turbulent flow for a given phase of a small-amplitude interfacial 
wave. In this approach, an eddy viscosity is formulated in terms of the typical scale of a 
turbulent eddy, which depends on the distance between the eddy itself and the air- water in- 
terface. Far from the interface, the turbulent eddies are advected quickly over an interfacial 
undulation, and have insufficient time to equilibrate, and so-called rapid-distortion theory 
is needed. Such an approach was developed by Townsend [15j. In the papers of Belcher 
and co-workers, however, the far-field region is simply modelled by the Rayleigh equation. 
Thus, the mixing-length is truncated: it is a simple function of the vertical co-ordinate close 
to the interface, and is set to zero far from the interface. We propose instead to follow the 
approach of Townsend [15j and lerley and Miles [2j: not only do we interpolate between 
the turbulent domains (a common factor between all these papers), but we also explicitly 
model the rapid-distortion region. Our model is based on the linearized closure model of 
the Reynolds- averaged equations pioneered by Launder et al. [L6\ and discussed in the book 
by Pope jnj. 

Numerical approaches: The work discussed so far uses asymptotic techniques for the de- 
termination of stability. Numerical methods provide more information on the stability, since 
they are valid over all parameter ranges, and since they readily enable the re-construction 
of the velocity and pressure fields from the solution of an eigenvalue problem. The paper 
of Boomkamp and Miesen provides yet more information, since they classify interfacial 
instabilities according to an energy budget: the contributions to the time-change in kinetic 
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energy are derived from the linearized Navier-Stokes equations, and several destabilizing 
factors are identified. Of interest in the present application are the critical-layer mechanism 
and the viscosity-contrast mechanism. In this approach [181 HH], and in similar work by 
Ozgen [201 EJ, the viscous liquid layer and gas layers are modelled co-equally; these added 
complications necessitate a numerical solution of the equation. The interfacial wave speed 
is then determined, along with the growth rate, as the solution of an eigenvalue problem. 
Turbulence enters the problem only though the logarithmic mean-fiow profile chosen, and 
the Reynolds stress terms are ignored. 

Modelling of the base flow: The papers mentioned so far use a model of the base fiow 
that captures its logarithmic shape. These base models contain a free parameter, however, 
since the interfacial friction velocity U^i is undetermined. In the approach we take, we 
use a closed model that fixes this parameter as a function of the Reynolds number. Our 
formalism is based on that introduced by Biberg [2^, although it differs in one key respect: 
we pay close attention to the near-interface modelling (and also to the modelling near the 
upper wall, if present), and remove the singularities that arise when the logarithm in the 
velocity is evaluated an upper wall or an interface. This is accomplished using a Van Driest 
interpolation ^23j , in which the velocity field smoothly transitions between the log layer 
and the viscous near-wall (interfacial) region. 

These are the ingredients of our model, which we develop in the following way. In Sec. [IT] 
we provide a detailed description of the rapid-distortion and near-equilibrium turbulence 
which we constitute in terms of the streamfunction. We introduce the Orr-Sommerfeld 
equation for the streamfunction, which couples to the turbulent stresses. In Sec. ITTT] we 



describe the numerical method used to solve the model equations. In Sec. [IV] we describe 
the effects of the turbulent modelling on a thin-film fiow (this problem was studied in the 
absence of perturbation turbulent stresses by Miesen and Boersma [19j). We find that the 
growth rate is not significantly affected by the turbulence modelling, although the structure 
of the fiow field is. Nevertheless, for deep-water waves (Sec. V), both the growth rate and 
the fiow structure are modified by the turbulence (the growth rate is enhanced). In this 
deep-water study, we compare our results with direct numerical simulation, and investigate 



a mechanism for generating thin-film type waves in a deep water setting. Finally, in Sec. (VT 
we present our conclusions. 



II. THEORETICAL FORMULATION 

In this section we provide a problem description and develop a turbulence model that 
takes account of rapid distortion and near-equilibrium regions in the turbulent stratified 
two-phase fiow. Our approach is inspired by the paper of Miles [2j. It takes into account 
the interactions between interfacial waves and the turbulence. The work here is, however, 
multi-faceted: there is a model of the fiat-interface turbulent state, together with a linear- 
stability analysis around this state, in which the effects of the turbulence are modelled in 
detail. The two-layer system we consider has the same structure as boundary-layer fiow, and 
is described schematically in Fig. [T] A rectangular co-ordinate system is used to model the 
fiow, with the plane z = coinciding with the undisturbed location of the interface. Since 
the interfacial waves that develop are two-dimensional, we restrict to a two-dimensional 
system, in which the bottom layer is a liquid of depth rf^,, which in this paper will either be 
a thin laminar layer, or a deep-water fiow (rf^ = co). The top layer is gaseous, turbulent 
and fully-developed. A shear stress is applied along the channel. The mean profile is thus 
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Figure 1: A schematic diagram showing the system of interest. There is turbulent shear flow in 
the gas layer, while the liquid layer either represents a thin laminar fllm, or deep-water flow, as in 
oceans. 

a uni-directional flow in the horizontal, x-direction. In the gas layer, near the gas-liquid 
interface and the gas- wall boundary, the flow is dominated by molecular viscosity, since here 
the viscous scale exceeds the characteristic length scale of the turbulence [TTl [24] . In the 
bulk of the gas region, the flow possesses a logarithmic proflle [17^ |24j. Finally, we assume 
that the gas-liquid interface is smooth. 



A. The perturbation equations 

Following standard practice dating back to Miles [6], we base the dynamical equations 
for the interfacial motion on the Reynolds-averaged Navier-Stokes (RANS) equations. The 
turbulent velocity is decomposed into averaged and fluctuating parts C7 = ([/, W) and 
{u\w') respectively. The averaged velocity depends on space and time through the RANS 
equations: 

V • t/ = 0, (lb) 
where (•) denotes the averaging process. The constants p and fi denote the density and 



viscosity respectively. In addition to viscous stress, Eq. (la) contains a stress term that 
arises due to the interaction between the turbulence and the mean flow. We use these 
equations to model a flat-interface or base state of the two-phase system shown in Fig. [T] 
Next, we introduce a small disturbance that shifts the flat interface at z = to z = ry, where 
1 77 1 <C 1, and this induces a change in the average velocity and pressure flelds, denoted as 
follows: 

([/, W, P) = ([/q {z) + 5u (x, z, t) , 5w (x, z, t) , Po z) + 5p (x, z, t)) , 



where we denote base-state quantities by a subscript zero. Since the flow is turbulent, and 
since the perturbations take the form of waves, they must satisfy the RANS equations for a 
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are the perturbation stresses due to the turbulence in the perturbed state, while the quan- 
tities with the zero-superscript are base-state stresses. Using the streamfunction represen- 
tation {bu^bw) = {(j)z^—(j)x)) and the normal-mode decomposition dx = ioi^ the perturbed 
RANS equations reduce to a single equation: 
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= p ip^ - a^Y (j) -\- piaT)5a + (D^ + a^) 5t, 



(3) 



where D = d/dz. 



B. Separation of domains 



To model the turbulent stresses in detail, we need to understand the features of the 
turbulence. Because the instabilities we observe arise due to conditions in the gas, we 
pay particular attention to the gas layer. If the liquid layer is deep, then the standard 
Orr-Sommerfeld equation gives an appropriate description of the flow there, provided an 
accurate base-state proflle is supplied jlHl [25|, while if the liquid resides in a thin layer, 
the Orr-Sommerfeld equation gives an exact description of the flow. Thus, in both cases, 
we specify the standard Orr-Sommerfeld equation in the liquid, without the additional 
stresses. The turbulence in the gas is characterised by two timescales: the eddy turnover 
timescale, and the advection timescale. Roughly speaking, the eddy turnover or turbulent 
timescale is the time required for a typical turbulent eddy to interact with the surrounding 
fluid and come into equilibrium (where the production equals dissipation). An estimate 
of this timescale, at a distance z from the interface, is Tt = Kz/U^i^ where k is the Von 
Karman constant, n = 0.41, and U^i is the friction velocity at the interface. The advection 
timescale is the time needed for the flow to advect an eddy over a wave crest. This flow 
distorts the turbulence and moves it away from equilibrium. An estimate of this timescale 
is Ta = [a\Uo (z) — c\]~^^ where a is the wavenumber and c is the complex the wave speed. 
These eddy-turnover and advection eflFects compete: near the interface, Tt is small compared 
with Ta, that is, the frequency of turbulent interactions is large compared with the frequency 
of advection events. This is the region of near-equilibrium where by deflnition eddy viscosity 
and one-equation turbulent closures are expected to be appropriate. Far away from the 
interface, Tt is large compared with the advection timescale (at least for the large, most 
energetic turbulent structures), or equivalently, the frequency of turbulent interactions is 
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small compared with the frequency of advection events. This is the region where rapid- 
distortion theory is expected to apply [121 HSl • obvious way, we call the region 
where eddy viscosity and one-equation formalisms are appropriate the near-interface region^ 
while the region where the rapid-distortion calculations work is called the far-field region. 
Crossover occurs at z = Zt^ for which Tt (zt) = Ts,{zt). For further discussion of this 
separation of domains, see ^141 126] . 



C. Turbulence modelling in the distinct domains 

In modelling the perturbation turbulent stresses (PTS), we are particularly interested in 
the anisotropy tensor 

where q = 2k = tr(rjj), and k is the turbulent kinetic energy. The exact equations for the 
fluctuating component of the velocity in the RANS equations are 

d \ , drik , dUi , du'i ^ , I dp' ^ 

In the far field, the mean strain produced by the wave is large. We therefore study the 
equations of motion in the limit of strong mean strain. In that case the terms involving 
gradients of fiuctuating quantities drop out: 

where 

p dxj dxi ' 

with p*^^^ being the so-called rapid pressure. By a straightforward manipulation of Eq. ([6]), 

^ + ^ (f/^r,, + T,,,) = p., + (7) 

in the limit of strong mean strain, where the tensor quantities have the following meaning: 
• The production of stress by the mean fiow: 



dxk dxk ' 



• The transport of stress by the rapid pressure: 



• The pressure rate-of-strain term: 



\ p \dxj dxi 
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For waves, the shear rate a{Uo — c) enters into the advection term in Eq. ([t]). In the 
rapid-distortion domain, this rate is large, and thus advective transport dominates over 
pressure-driven transport. We therefore omit the pressure transport term in ([T]), which now 
simphfies further: 

| + t/-vV..=^..+^?, (8) 



Equation (|8|) contains only one term that is not available in closed form: the rapid pressure- 



rate of strain tensor This term has been modelled accurately by Launder, Reece and 

Rodi [16j (see also [17]). They use the following form: 

1^ = -Ci {v., - lV,,S,,) , (9) 

where Ci = |. This model also has the property of linearity in the tensor r^^, which is 
desirable in any rapid-distortion theory [17j. The advection equation ([s]) now becomes 



-^ + U- Vxi, = (1 - Ci) Vi, + \CiVkAr (10) 
Using the definition ^ and Eq. (10), the n^j's satisfy the relation 

^ + C/ . Vn., = - (1 - Q) n.,g - (1 - Q) n,,g - 2n,,g {\CA, - n.,) . 

We decompose the tensor n^j into a flat-interface component and a perturbation component: 
Wij = wf^ + 5wij^ where wf^ is the equilibrium anisotropy tensor, and is thus constant in 
time. Hence, 



ni2 = - — + wf^Sq) , 
nil = nSi^ [5a^ + nji^^g ) , 

(7n V / 



where {z) represents twice the turbulent kinetic energy associated with the base state 
and 5q a small-amplitude perturbation about this state. We therefore have perturbation 
equations for 5wi2 and 8w := ^rin — 5r\22'' 

d \ d d d 

— 7 + C7 • V 5wi2 ai—5u + a2—5w + a^—5u + a4^5wi2 + a^5n, 
at ) oz ox ox 

|- + (7 • V ) (5n = (3i^5u + l32^Sw + f3^^5u + [3^5x\i2 + /^s^n. (11) 
ot I oz ox ox 
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Given the base-state shear rate S = dUo/dz^ the coefficients in Eqs. (11) are defined as 
foUows: 
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Far from the interface, the effects of rapid distortion due to the base fiow are neghgible, 
since the shear rate S (z) is smaU there; thus 0^4, 0^5, /34, and f3^ are smaU and are neglected 
in this formahsm. In Sees. 



IV 



and 



V 



we estimate the n^^^'s by constant values obtained 



from the properties of log layer in turbulent shear fiow. Then, Eq. (11), in normal-mode 
form, reduces to a set of equations that are algebraic in 5ni2 and 5n: 



ia {Uq — c) 5ni2 
ia {Uq — c) 5n 



ai—du + a2—ow + as—du, 
oz ox ox 

d d d 

oz ox ox 



(12) 



The right-hand side of both these equations can easily be recast in terms of the stream- 
function, using {d/dz)5u = D'^cj) etc. However, to highlight the effect of the distortion of 



the Reynolds stresses by the perturbation velocity, we leave Eq. (12) in its present form 
throughout the paper. 

The derivation of a near-field model is less involved. Near the interface, the turbulence 
tends quickly to equilibrium, and the n^j's take their constant value. Thus, 

5n,j = 0. (13) 



Next, we develop a strategy to join up these domains and the models (11) and (13) into a 



single formalism. This requires a detailed picture of the turbulent kinetic energy. 



D. Turbulent kinetic energy 

The governing equation for the perturbation turbulent kinetic energy Sk = ^{u^'^+w^'^) — ko 
can be written as: 

^^5k + Uo • V5k + 5U 'Vk = -V ' {5 J) + 5V - fe. 
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where the terms are explained in detail in what follows. Using the properties of the basic 
state, this equation reduces to 

—5k + Uo^Sk + 5w-^ = -V • (5 J) + 5V- 5e. 
ot ox dz 

We now turn to the problem of modelling the terms on the right-hand side of this equation. 

• The energy flux: In Townsend [151 EIJ and lerley and Miles [2j, the energy flux 5 J is 
ignored, since turbulent transport is negligible in the bulk gas flow. It is important, 
however, near the interface, where molecular viscosity dominates. We therefore take 
5 J = -uVSk. 

• Dissipation: In previous works, [21 [151 [27], the dissipation function 6s in the bulk gas 
region was constructed as a cumbersome function of ko and the mixing length. In this 
work, we make use of the simple model 

6e = —6k^ 
X 

where x is a timescale constructed from dimensional analysis: 

V 

^ *i 

We have verifled that the choice of form for this relation makes little diflFerence to the 
results of the stability calculation, and we therefore settle on the linear form, which 
has the added advantage that does not diverge at z = 0. 

• Production: The production of turbulent kinetic energy is half the trace of Vij^ and 
thus the perturbed production rate has the form 



5V = -5u,^^-rf^^5U,, 



'•^ dxj "-^ dxj 



which can be re-expressed as 



where St = —{u'w') — r^^^ is the perturbed turbulent shear stress, and t^^\ (j^x \ and 
(jf^ represent the turbulent stresses in the flat-interface state. Thus, we obtain the 
following kinetic-energy equation: 

ia {Uq — c) 6k + x~^^k 

= yCD^- a^) 5k - Sru^ - rg^ (D^ + a') 4> - mr^Dc/) + ia^4>, 
^ ^ dz ^ ^ dz 



where r^^^ = — ^ai^^ — ai^^^ and vf^2 = 



-T 



(0) 
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E. Interpolation between domains 



The algebraic model for the tensor 5n^j is given by Eq. (12) in the far field, and by Eq. (13) 



in the near field. We can interpolate between these domains by using a hybrid stress model: 

f d d d \ 

[ia {Uo - c)I{z) + (1 -I{z))] 5ni2 =1{z) ( ^i^^'^ + + ^^~q^^^ j ' 

/ C\ o o 

[la {Uo -c)l{z) + {l-l {z))] 5n=l {z) (f^iQ-Ju + [^^^Sw + [3^—5u 

where Z {z) is an interpolating function that is zero at z = and asymptotes to X = 1, with 
a characteristic lengthscale Zf 

We now have all the components to assemble the model, which we present below in non- 
dimensional form. We non-dimensionalize on as- yet undefined length and velocity scales 
h and V respectively, and the gas viscosity, which gives rise to the gas Reynolds number 
Re = Vh/i^G- In the liquid, a single equation describes the fiow, which in non-dimensional 
form is simply 



{Ul - c) (D^ - a') 



m 
Re 



2\2 



(D2 - a') 



z<0, 



(14a) 



where the subscripts 'L' and 'G' denote evaluation in the gas and liquid phases respectively: 
Ul and Ug are the base-state velocities in the liquid and gas phases, and r = Pl/pg and 
^ = Pl/pg are the density and viscosity ratios. In the turbulent gas, we must solve the 
following coupled set of equations: 



(Ug - c) 



a^) (j)G 



dz^ ' 



1 

'Re 



(D^ - o?y (j)G - iaBSr - (D^ + a^) Sr^, 

(14b) 



ia {Ug - c) + 



Rel 
Re 



Sk 



^ (D2-a2)5fc_5rl2^-rS?(D2 + a2^^^_,«r(0)D(/> + ^a^(/>G, (14c) 
^ ' dz ^ dz 



Re 



where i?e* = Unh/fc, the Reynolds number based on the interfacial friction velocity 
and where 

(0) 

[ia {Ug -c)X{z) + {1-X {z))\ Sr,^ - [ia (Ug -c)X{z) + {l-X {z))] 5q 

f d d d \ 

= qo {z)X{z) iai—Su + a2-^5w + 03— (5« j , (14d) 



[ia{UG-c)X{z) + {l-X{z))] Sr 



r(o) 
Qo 



[ia{UG-c)Xiz) + il-X{z))]Sq 



qo{z)X{z)[(3i—8u + (32^5w + (3^—8u\ . (14e) 
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At the waU z = —di/h, we have the no-shp conditions on the streamfunction (j)L, 

4>L i-dL/h) = D0i i-dL/h) = 0. 

The equations are also matched across the undisturbed interface z = 0, where we have the 
following conditions: 

(f)L = 4>G, (15a) 

m (D^ + a^) = (D^ + a^) + ReSro, (15c) 

m (D^(I)L - 3a^D(l)L) + iarRe (c - Ul) ^(t)L + iarRe—^(t)L — (Fr + a^S) (t)L 

^ ^ dz c — Ul 

= (D^(t)G - 3a^D(t)G) + iaRe (c - Ul) B(f)G + iaRe^^(t)G + ReDSr + iaReSa, (15d) 
^ ^ dz 

where we have introduced the inverse Froude and inverse Weber numbers, respectively 

Here g and a are the gravity and surface tension constants. For clarity, we shall refer to 
the inverse Froude number as the gravity number, being the ratio of gravitational to inertial 
forces, likewise we refer to the inverse Weber number as the surface-tension number. 

The turbulent kinetic-energy equation is second-order, and therefore requires only two 
boundary conditions. We apply the conditions ko + 5k = on z = r] and z = dG/h^ 
where rfc/^ is the non-dimensional vertical extent of the gas layer. For boundary-layer 
flow, rf^ = oc, while for channel flow rf^ is flnite. Upon linearization these kinetic energy 
conditions are 

6k = at z = 0, and at z = rf^/Zi. (17) 

Furthermore, if we take 

dX(0) ^ 
1(0) = ^=0, 

then, by inspection of Eqs. (14d) and ( |14e ), Sri2 = {fu/l) = at z = 0, and similarly 
for Sr and DSr. (The justification for this choice of X is given at the end of Sec. IV A ) 
Consequently, Eqs. (fT5|) reduce to the standard interfacial equations: 



<j)L = (j)G, (18a) 

m (d2 + a^) = (D^ + a^) (Pc, (18c) 

dUr ICXT R,c 

m (D^4)L - 3a'^D4>L) + iarRe (c - Ul) D4>l + iarRe— —4>l — (Fr + a^S'cap) 4>l 

^ ' dz c — Ul 

= (dVg - Sa^Dcpo) + iaRe (c - Ul) B(j)G + iaRe'^cpo- (18d) 
^ d^ 
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Figure 2: A schematic description of the left-hand side of Eq. (20). The shaded blocks represent 



the linear interaction between the turbulence variables. The upper right-hand block enclosed by 
the thick dotted lines is the matrix for the eigenvalue problem in the absence of the turbulent stress 
and energy variables. 



Finally, at the top of the gas domain z 
tion velocities should vanish: 



dc/h^ we have the requirement that the perturba- 



{dclh) = D(t)G {dc/h) = 0. 

Thus, we obtain a total of 10 conditions, which provides sufficient information to close the 
system of equations {(l)L^(pG^^k^5vi2^Sv). No conditions are required on Sv^ and 5ri2, since 
these variables appear only in an algebraic way in Eqs. (14d) and (14e). We now turn to 
the numerical solution of this system of equations 



III. NUMERICAL METHOD 

In the liquid and turbulent gas domains, we propose an approximate solution to the 



Orr-Sommerfeld (OS) equations (14) in terms of Chebyshev polynomials, such that 



Nl Ng 

(z) ^ ^ a^Tn (tjl) , (z) ^ ^ bnTr, {tjg) , (19a) 

n=0 n=0 

Ng-2 Ng-4: Ng-4: 

Sk=Y, CnTn iVG,i) , 5ri2 = J] 9^T^ (tig) , Sr=Y^ inTn (vg) , (19b) 

n=0 n=0 n=0 



13 



where we have mapped the domains of the various layers into the interval [—1, 1]. For our 
two-layer system, with z G [— ^l, 0] for the liquid and z G [0, rf^] for the gas turbulent layer, 
the appropriate coordinate transformations are the following: 

Thus, we assume that the boundary layer is of finite depth. This setup has the advantage 
that the gas domain is finite in extent, and thus there are no virtual interfaces in the problem. 
In our experience, these virtual interfaces can interfere with the stability properties of the 
system, and are therefore undesirable. Although there are differences between the stability 
behaviour of the finite- and infinite-domain frameworks, these are apparent only at small a 
(large- wavelength disturbances) . We shall also show that this is far from the most dangerous 
mode, and these differences are therefore unimportant. The derivatives of our approximate 
solution are obtained through the expressions 

'^-hA'^^n^.o), (190) 

n=0 ^ ^ n=0 ^ ^ 

where the derivatives of the Chebyshev polynomial (ry) are obtained through the appro- 
priate recurrence relation. 



To reduce the problem (14) to a matrix form, we substitute the approximate solution (19) 



into the OS equations (14) and evaluate the resulting identity at the collocation points 



= cos 2 j ' V=^^ - 3, 

m = cos (^ ^^^_ 2 ^ ' P=^^ - 3. 

In the two-layer application that we have in mind, this gives 

{Nl - 3) + 4 (TVg,! -3) = Nl + 4iVc,i - 15 

equations in (A^^ + 1) + {Ng,i + 1) + ( Ag,i - 1) + 2 {Ng,i -S) = Nl + 4:Ng,i - 5 unknowns. 
We must therefore find an additional 10 equations to close the system. These are provided 
by the 10 boundary and interfacial conditions. Thus, we are reduced to solving for the 
wave speed c, given + 4:Ng^i — 3 linear equations in as many unknowns. In symbols, we 
therefore have a generalized eigenvalue problem 

Lx XMx, A -iac, (20) 

where positivity of the real part of A (Ar > 0) indicates instability. A schematic description 
of the left-hand side of this system is shown in Fig. [2| The solution of (20) is best found 



using a linear algebra package (as in [28]). Then, a built-in eigenvalue solver automatically 
balances the matrices L and M, an important issue here since these matrices can be very 
badly conditioned, owing to the high-order derivatives of the Chebyshev polynomials that 
appear in the expansion of the streamfunction (as in the work by Boomkamp et ai [29]). 
Furthermore, a specialized eigenvalue solver takes into account the appearance of zero rows 
(hence infinite eigenvalues) in the matrix M, and thus gives an accurate answer. Following 
standard practice, we modify the number of collocation points until convergence is achieved 
(typically, = 200, Nq = 150 for convergence to four significant figures). 



14 



IV. RESULTS FOR A THIN FILM 

In this section, we investigate the effects of turbulence on the interfacial stabihty of a 
thin laminar liquid layer exposed to fully-developed turbulent shear flow in an overlying gas 
layer. The parameter regime we study involves i?e, Fr, and the ratio 




liquid-fllm thickness 
gas-layer thickness 



and is chosen to mimic the conditions in pipe flow in oil-gas transport problems [19j. The 
liquid Reynolds number is such that the liquid flow is laminar. We use a base state that 
mimics flow in a boundary layer. In this case the flow is conflned by a flat plat at z = rf^, a 
large distance from the interface. This plate moves at velocity Uq relative to the interface. 
It is thus appropriate to re-scale on the velocity [/q, and on the gas-layer depth do (thus 
h = do and V = Uq^ in the notation of Sec. [ll]). This setup has the advantage that the gas 
domain is flnite in extent, while having little effect on the results (see Fig. [5] for a comparison 
between flnite- and inflnite-domain growth rates). Using this framework, we constitute the 
base-state velocity. 



A. Base-state determination 

The gas layer: In the gas, we have the following momentum balance for the base state: 

I^G-^ + ttss = (21a) 

where ri is the interfacial shear stress, and ttss is the turbulent shear stress —pg{u'w'). To 
close this term, we make use of the interpolation function 

G (5) = 5 (1-5), 5=-^, (21b) 

do 

where z is the non-dimensional vertical co-ordinate (the tilde decoration has been introduced 
to designate dimensionless quantities). Now the turbulent shear stress has a characteristic 



Taylor expansion near the walls, and for the interpolation function (21b) to agree with this 
expansion, we must modify Eq. \21h\ near the upper wall and the interface, such that 

TTSS = KpcdcU^iG (5) ^ (5) ^ (1 - 5) ^ (5) = 1 - e-^'"/^, (21c) 

where n and A (which depends on Re) are model parameters flxed below with reference to 



Eq. (22), and = ^jT\j pc is the interfacial friction velocity. Hence, 



Ug = Tide 







where we have set the frame of reference to move with the interface. We non-dimensionalize 
on the plate velocity [/q, on the gas height rf^, and on the gas viscosity /xg, which gives a 
non-dimensional velocity 
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where Re^ = pcU^idc/ I^^g^ in contrast to the Reynolds number set by the non- 
dimensionahzation, Re = PgUqcIg/ /^^g- The Reynolds number Re^ is determined as the 
root of the equation 

M f = 1 (21e) 

Re Jo 1 + KRe^G (s) xP{s)xP{l-s) ' ^ ^ 

The liquid layer: The momentum balance in the laminar liquid is simply 

/..^ = r. (21f) 
Integrating and non-dimensionalizing, this is 



The turbulence-related variables: From Eq. (21c), we obtain the following non- 
dimensional form for the gas Reynolds stress: 

_ ^r^^ _ Rel kRc^G (z) {z) {1 - z) 
^^^^ ~ pgU^ ~ Re^l + ^Re,G (z) ^Jj (z) ^Jj {1 - z)' ^ ' 

Moreover, we constitute the turbulent kinetic energy as 

k 1 Rp'^ 

where C is another constant. We take C = 0.55, which is the value appropriate for the 



logarithmic region of the mean velocity in a boundary layer. The form (21i) takes into 
account the constancy of the turbulent kinetic energy in the core, while damping the energy 
to zero at the interface and at the wall. Well into the bulk of the gas flow, this model gives 
the ratio 

T 



(0) 

2k ' 2 



— n^^ — — ^ 



To model the quantity r^^^ = r[^^ — = {u''^) — {w'^) (where this average is over base-state 
variables), we make the approximation 

ri? = aii/c, r^22 = ^22^, (21j) 

where an and are constants, such that an + a22 = 2. For linear algebraic stress models, 

the convention is to take an = a22 = 1, which is physically incorrect: the turbulent kinetic 

energy is not equipartitioned in the streamwise and normal directions. One flx is to use a 

non-linear algebraic stress model [30], which gives better predictions of the partitioning of 

energy. For simplicity, however, we use typical values for the partition constants from DNS 

data in boundary-layer flow [31j. Note flnally that this linear relationship, while true in the 

bulk flow, breaks down near the interface and near the wall. We have, however, conducted 

tests to verify the sensitivity of the stability analysis to this detail and have found that 

3 



changing the form of the law (21j) has little effect on the results. Thus, we take an = ^ 
throughout the flow, and hence, a22 = |. 

rW = fc, nff = |, ng) = i (21k) 
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Figure 3: The base-flow variables for Re = 8000, m = 55, r = 1000, and 5 = 1/50. Subfigure (a) 
shows the velocity profile in the liquid and the gas, with distances normalized on the liquid-film 
thickness; subfigure (b) shows the development of a logarithmic profile near the centreline of the 
gas domain. Subfigure (c) shows the spatial dependence of the turbulent shear stress and kinetic 
energy. 

Close to the interface 5 = 0, we have the following asymptotic expressions for the model: 

~i+n ] f i^^Re^\ fj ^ kRc^ ^^^2 



The asymptotic behaviour of these model variables can be made to agree with the true 
behaviour of near- wall turbulence. Although the system we study possesses both a wall 
and an interface, when the density contrast between the liquid and the gas is large, DNS 
results suggest [32j that a comparison between wall- and interfacial-turbulence is justified. 
The fluctuating velocities for wall turbulence have the following Taylor expansions around 
the wall location z = 0: 

u' = ai (x, t) + hi (x, t) z + Ci (x, t) + w' = a2 (x, t) + 62 t) z + C2 (x, t) + .... 

Using the no-slip conditions, we obtain ai = a2 = 0. Since is zero along the wall z = 0, 
{du' / dx) 0, hence {dw' /dz)^^^ = 62 = 0. Thus, near the wall, 

k = \pG{u'' + w'')r.\pa{hl)z\ 

T = -Pg{u'w') - -pG{hrC2)z\ (22) 

Note that this result provides for a no-flux condition on the turbulent kinetic energy k at 
the wall, dk/dz = on z = 0. Choosing n = 2 in our model forces agreement between the 
interface model and the near- wall asymptotics. The results of the model are shown in Fig. [3j 

Computation of Zti In Sec. [n| we introduced the turbulent and advective timescales. 
These timescales vary in the normal direction. Close to the interface, the frequency of eddy 
turnovers is T^~^, which is large compared to the advection frequency Ta. Far from the 
interface, the magnitude of these frequencies is reversed. Crossover occurs at Zt, which is 
determined as the root of the equation 
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Typically, Cr ^ and negligible error is incurred by replacing c by in Eq. (23). Making 
this replacement, we can readily find a lower bound for Zt\ 

Re^a~^ = nzt Re \ U {zt) — Cr| , 

< KZt Re U {zt) , 

< K.Relzl, 



hence 



(24) 



When the wave speed is small, the bound (24) sharpens. Moreover, we have verified that 
the propagation speed is not aflFected by the turbulence modelling (e.g.. Fig. 11 (b)), and 
thus the computation of Zt can be carried out using the model without the perturbation 




Figure 4: The crossover height zt as a function of wavenumber for Re — 8000. The wavelength 
is shown for comparison. 

turbulent stresses (PTS). In Fig. [i] we have done precisely this: the scale Zt is computed 
and the bound ([24]) is found to be sharp. The scale which determines the spatial 
extent of the streamfunction, is shown for comparison. For all but the largest a-values 
(shortest wavelengths), the crossover height Zt lies in a region where the streamfunction is 
non-negligible. Thus, we expect the streamfunction to 'feel' the effects of rapid distortion. 
Note finally that the computation of the scale Zt enables us to develop an interpolation 
function X(z), specifically, we take 

l{z) = l-e-("/"^)', (25) 

in the gas. The exponent is taken to be 2 to guarantee that the turbulent eflFects vanish 
strongly at the interface: since the fiat-interface state has this property, it is reasonable to 
assign the same property to the perturbed state. 



B. Linear-stability analysis 

We carry out a stability analysis around the base state just constituted for Re = 8000, 
Fr = 10~^, and S = 10~^. The density and viscosity ratios are chosen such that the stability 
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Figure 5: (a) The growth rate in the thin-film case, with Re — 8000 (sohd hne). The wavenumber 
is normahzed by the film thickness d^,, rather than the gas-domain thickness d^, and dc/dL — 50. 
The difference between the PTS wave and the non-PTS wave is so small as to be negligible, and 
is not shown here; (b) The wave speed in the thin- film case (solid line). A comparison with a 
boundary-layer base state with no upper plate is shown (dashed-dotted lines). 



analysis models an air-water system under standard conditions m = 55, r = 1000. These 
Reynolds, inverse Froude, and inverse Weber numbers are chosen to correspond to typical 
air- water conditions for a liquid-film thickness d^ = ^^g/SO. We examine the growth rate 
of the disturbance with and without the effects of the PTS, and the results are shown in 
Fig. [5| where we also plot the wave speed. The inclusion of the PTS has little effect on 
the growth rate in this case, except at long wavelengths. The trend in the dispersion curve 
for the PTS wave shows a tiny downward shift in the growth rate relative to the non-PTS 
wave. The maximum growth rate occurs dX a 60. Since we have non-dimensionalized 
on the gas-layer height, this corresponds to a wavelength l/dL ^ 2Ti5/a = 27r x (50/60), 
that is, a wavelength greater than, but comparable to the liquid film thickness. Finally, we 
note that our finding that the maximum growth rate occurs at a wavelength comparable 
to the liquid- film thickness is similar to results obtained elsewhere for thin liquid films. 
The U-shaped trough in the wave-speed diagram is also characteristic of thin-film flow: we 
have observed it by reproducing the work of Miesen and Boersma fT9] , and extracting the 
wave speed from the stability analysis, in addition to the growth rate. Our inclusion of an 
upper plate has little or no effect on the results: this can be seen by comparing the data 
in Fig. [5] with a true boundary layer, at the same parameter values (in particular, at the 
same liquid Reynolds number). This is shown by the broken-line curves in Fig. [sj where the 
boundary-layer stability analysis is conducted according to the framework of Miesen and 
Boersma p!9]. 

Figure [5] implies that the PTS have little effect on the growth rate at the Reynolds number 
considered {Re = 8000). This can be explained by studying the form of the rapid-distortion 
equations ( |14d[ ) and ( |14e[ ) : here the stress terms are damped not only by the decay function 
X(z/zt), but also by the equilibrium turbulent kinetic energy function q{z)^ which is zero 
at the interface. For now, let us take it on trust that the instability in Fig. [5] is due to the 
viscosity-contrast mechanism (we verify this in Tab. IT]). Then, the instability is governed 
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by conditions at the interface and in the viscous sublayer, z < Re~^. This is precisely the 
region in which q (z) is damped to zero. Thus, in the region z < Re~^ where the growth rate 
is determined, rapid distortion plays no role, and therefore, does not affect the growth rate. 
A forteriori^ the PTS play no role in determining the growth rate for this particular class 
of instability. Nevertheless, there exists a further aspect of the problem where the PTS do 
alter the features of the disturbance significantly, namely the shape of the streamfunction 
in a zone far from the interface. Indeed, we shall find that the far-field streamfunction 
has an oscillatory structure that propagates far into the gas core. The existence of these 
oscillations casts doubt on whether the structure and statistics of the turbulence near non- 
linear interfacial waves is similar to that for a wall region. Ultimately, this must be confirmed 
by accurate numerical simulation. 

We investigate the spatial structure of the perturbation velocity for wavenumbers a = 15 
and a = 60] the results are shown in Figs. [6] and [7| The eflFect of the rapid distortion on the 
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Figure 6: Flow field structure for a = 15. Subfigures (a) and (b) show the streamwise velocity field 
and a profile at x = 0; subfigures (c) and (d) show the normal velocity field and a profile at x = 0; 
subfigures (e) and (f) show the pre-averaged version of the wave Reynolds stress, namely the field 
uw, and its profile at x = 0. In each case, we have normalized the velocities by max \u\. 
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Figure 7: Flow field structure for a — 60. Subfigures (a) and (b) show the streamwise velocity field 
and a profile at x = 0; subfigures (c) and (d) show the normal velocity field and a profile at x = 0; 
subfigures (e) and (f) show the pre-averaged version of the wave Reynolds stress, namely the field 
uw^ and its profile at x = 0. In each case, we have normalized the velocities by max \ u\. 



flow is visible at a = 15, and almost negligible at a = 60. For a = 15, the streamfunction 
extends into the rapid-distortion domain, and the wave and the turbulence interact. This 
gives rise to an oscillatory structure in the velocity field, in the normal direction, as shown 
in Figs. [6] (b) and (d). This structure also manifests itself as streaks, shown in Fig. [6] (a), 
which are inclined at an acute angle relative to the interface. These oscillatory structures 
are all but invisible in the a = 60 case (Fig. J7|. 

We plot the pressure distribution in Fig. |8j The minimum pressure is almost exactly out 
of phase with the interfacial variation at a = 15, while at a = 60 the pressure minimum is 
shifted slightly downstream of the free-surface maximum. The co-incidence of the pressure 
minimum and the free-surface maximum is explained by Bernoulli's principle, since the gas 
viscosity is small. The in-phase component of the pressure at a = 60 is a consequence of 
the so-called quasi-separated sheltering O [29] , wherein the viscous sublayer creates a wake 
downstream of the crest, which reduces the pressure fluctuation there. 
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Figure 8: The pressure distribution for a — 15, (a), and a — 60, (b), normalized by max|i/|, the 
maximum streamwise perturbation velocity. 

The other phase relationships, namely those between the velocity at the interface and the 
interfacial disturbance can also be understood quite readily. The normal component of the 
velocity possesses successive extrema (maxima or minima) . The first set of extrema is located 
at the surface z = 0, and can be understood by examination of the kinematic condition, 
which in the moving frame of reference is simply drj/dt = w dit the interface. Upstream of 
a wave crest that is propagating from left to right (cr > 0), the interfacial height decreases 
in time, which implies that the velocity w is negative there, while downstream of the wave 
crest, the vertical velocity component is positive. Moreover, since w = 
u = ^ |'g2a(x-ct) / (iz)] ^ there is a 7r/2 phase difference between u and and thus the 
structure of the streamwise velocity can be understood simply as a phase shift relative to 
the normal velocity, in order to satisfy the continuity condition. 

Finally, we make use of an energy-decomposition to pinpoint the source of the instability. 
This decomposition or budget is obtained from the RANS equations, and was introduced 



a = 15 


KINl 


KINg 


REYl 


REYg 


DISSl 


DISSg 


TURB 


NOR 


TAN 


PTS 


0.97 


0.03 


0.02 


9.21 


-0.20 


1.35 





0.000 


1.000 


No PTS 


0.005 


0.012 


-0.001 


-0.185 


-0.019 


-0.738 


-0.039 


0.000 


1.000 




a = 60 


KINl 


KINg 


REYl 


REYg 


DISSl 


DISSg 


TURB 


NOR 


TAN 


No PTS 


0.008 


0.011 


-0.002 


-0.089 


-0.024 


-0.863 





-0.004 


1.000 


PTS 


0.007 


0.011 


-0.002 


-0.103 


-0.023 


-0.832 


-0.018 


-0.004 


1.000 



Table I: Energy budget for the thin-film flow, normalized on the tangential term. 
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by Boomkamp and Miesen [TH]: 
d 



(26a) 



—Sn 



—Sn 



-5rn + Sr22 -Sri2 
-5ri2 



(26b) 
(26c) 

which we multiply by the velocity Suj and integrate over space. We obtain the following 
balance equation: 



V • -Uj = 0, 



KIN,- = REYj + Y DISS,- + ^ TURB,- + INT, 

j=L,G j=L,G j=L,G j=L,G 



(26d) 



where 



KIN, = \jj dx j dzr,5ul 



REY, = -r, 



j dx j dz Su^Swi'^^ 



DISS,- = I dx I dz 



dz 

2 ( ^Su,) + (^Suj + ^Sw^ + 2 (-^Swj 



TURB,- = Sj,G<r 



dz 



dx 

— ( 
dx 



d f d d 

oz ox 



(26e) 
(26f) 

(26g) 

(26h) 



Additionally, 



INT = J dx [5ul5Jl,zx + SwlJl,zz]^=o - J dx [5ugSJg,zx + ujgSJg,; 

which is decomposed into normal and tangential contributions, 

INT = NOR + TAN, 

where 

NOR = J dx [5wlSJl,zz - SwgSJg,zz]^=o , 

and 

TAN = J dx [SulSJl,zx - SugSJg,zx]^=o • 

The results of this study are given in Tab. [T| We find in each case that the tangential term 
is the source of the instability. This can be re-written as 

TAN = [/^(O) (m- 1) / 5T^,{x,0)r]{x)dx, 

Jo 

which is positive when m > 1, and when the viscous shear stress 5Txz 0) is less than 7r/2 
out-of-phase with the interfacial variation. Thus, the instability is driven by the viscosity 
contrast m > 1. 
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C. A simplified calculation to elucidate the rapid-distortion effects 



To understand the effects of rapid distortion on the structure of the streamfunction, rather 
than on the stabihty of the system, we resort to a highly simpUfied model, where analytical 
progress is possible. This model possesses as its solution the oscillatory structures visible in 
Figs. |6] and [7} We perform a linear-stability analysis around the base state 



U 




-l<z<0, 

0<z<b, 

z>b. 



(27) 



We make use of the following Orr-Sommerfeld and PTS equations: 



iar {mz - c) (D^ - a^) 4>l = ^ (D^ - a^) 



Re 



-1< z<0, 



ia {z - c) (D^ 



2\2 



Re 



0<z<b, 



ia {b - c) - a^) (1>G = ^ (D^ - a^) + (D^ + a^) r, z>b, 



(28a) 

(28b) 

(28c) 
(28d) 

where now Qoo, the constant turbulent kinetic energy in the far field, can be thought of as 
parametrizing the PTS. Given a set of boundary and interfacial conditions, it is possible to 
obtain a closed-form solution to this set of equations, in terms of exponentials and integrals 
of Airy functions [33', "M]. Our goal here, however, is simply to elucidate the oscillatory 
nature of the solution found in the full analysis. To that end, we take a characteristic value 
of the wave speed c, and obtain a solution to the system (28) in the far field z > b. There, 
the solution is (f>G = e''^, where <; solves a fourth-order polynomial equation: 

(1 + iaiReq) {<;/a)'^ - a^Reg {<;/af + [i {ai + ^2) Reg - 2 - iReb] is/af 

— asReq {<^/a) + iReqa2 + 1 + iRef, = 0, (29) 



where 



_ Req^ 
a{h-cy 



Re {b - c) 



a 



When 6 |c|, we have the condition li^e^l ^ \Req\^ and we can treat \Req\ as a smaU ex- 



pansion parameter. The lowest-order solution of Eq. (29) is then = a'^ or a'^+iaRe {b — c), 
and the rapid-distortion effects appear at first order. In fact, 



^1 -a 
^2 = -a^/^Re^ 



1 + 



goo aiRet 
{b-cf 2 



|i?eJ > 1. 



(30) 



Note, however, that the rapid-distortion eflFects appear non-perturbatively in ^2 for suffi- 
ciently large values of li^e^l^oo- The exact condition is o^i^oo Reb/{b-cy = 0(1), which 



equates to 



a < 



aiQooRe 



(31) 
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Using typical values for the thin-film waves {Re = 10^, ^oo = 10~^, \b—c\ = 10~^, ai = 10~^), 
we expect rapid distortion to affect the streamfunction at lowest order when a < 10. This 
analysis is verified by the streamfunction plots in Fig. [9| using = 10~^, where we take 
b = 0.25, and Re = 8000. The coefficients o^i, 0^2, and are taken from Sec. |IIC[ while 
the c-values are taken from Fig. [5] at a = 10 and a = 60. In Fig. M (a) we compare the 
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Figure 9: Far-field streamfunction for the simplified model, Re = 8000, b = 0.25, and c taken from 
Fig. [5j The wavenumber is o: = 10 in (a) and a = 60 in (b). At longer wavelengths, the rapid 
distortion induces an oscillation in the streamfunction, seen in the comparison between the ^oo — 
and ^oo = 10~^ curves in (a). At shorter wavelengths, this effect vanishes (see (b)). This agrees 
qualitatively with the streamfunction pattern in the full model. 



far-ffeld streamfunction associated with the mode (^2 with and without the PTS at a = 10. 
There is a clear difference between the two cases, while the difference between the two modes 
associated with is negligible (the difference is small and the plot is not shown). This is 
consistent with the analysis in Eqs. (30) and (31). The consistency extends to the a = 60 



case (Fig. [9] (b)), where there is little difference between the ^1 and ^2 streamfunctions. 
The enhanced oscillations in the streamfunction visible in Fig. [9] (a) are similar to those 
obtained by Zaki and Saha [35j in their study of the interaction between disturbances in 
the free stream and the boundary layer. There, the authors used the continuous-spectrum 
Orr-Sommerfeld equation as a model, and the qualitative agreement between our findings 
and theirs further vindicates our results. 

Note finally that there is another way for the rapid distortion to make itself felt at 

:^r, the speed b 

> \Reb\ ~ 



c is small, and 



zeroth order in the polynomial equation (29). When b 

as an expansion parameter, we obtain the zeroth-order 



\Ren 



By treating |i?e5 



polynomial to solve: 



In this case, the oscillatory term in the streamfunction (j) = is due mostly to the 
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rapid distortion, and the contribution from the osciUation coming from the base flow is a 
perturbation. This is made manifest when we write down the closed-form solution of the 
equation that exists when = 0: 




We expect this last mechanism to be important for fast waves, that is, for a case in which 
the equation b — (a) = has a solution. We therefore turn our attention to a situation in 
which such fast waves occur. 



V. RESULTS FOR DEEP-WATER WAVES 



In this section we investigate the effects of turbulence on the stability of an interface 
separating a deep body of liquid from a gas layer that acts on the interface by a fully- 
developed turbulent shear flow. 



A. Base-state determination 



As before, we use a base state that mimics flow in a boundary layer, and thus, the flow is 
conflned by a flat plate at z = rf^, a large distance from the interface. This plate moves at 
velocity Uq relative to the interface. Using this framework, the basic velocity is constituted 
as before: the non-dimensional velocity is given by Eq. (21d), where Re^ = pcU^idc/ I^^g^ and 



Re = pcUo/dG' The eddy- viscosity and wall functions G and ^0 have their usual meaning, 
given by Eqs. ( |21b ) and (21c) respectively. The sole difference between the proflle in Sec. IV 
and that used here is in the liquid, where we make use of the following deep-water profile [25] : 



UL = a {e'' - 1) 



(32a) 



where a and b are constants to be determined. There is in fact, only one constant to 
determine, since the continuity of tangential stress requires that 



m- 



dz 



dz 



0. 



(32b) 



Hence, mab = Rel/Re^ and there remains a single free constant a whose value is fixed in 
Tab [11} Thus, the liquid velocity has the reduced form 



Ut= a \ e^«^e — 1 



Computation of Zt: In carrying out the stability analysis, we have verified that the real 
part of the wave speed is affected only slightly by the turbulence modelling (less than %1), 
and thus the computation of Zt can be carried out using the model without the perturba- 
tion turbulent stresses (PTS). The crossover value Zt where the turbulent and advection 
timescales are equal is thus determined by Eq. (23); the dependence of Zt on a is shown 
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Figure 10: The crossover height as a function of wavenumber for Re 
m = 55. The wavelength is shown for comparison. 



10^ r = 1000, and 



in Fig. 10, for Re = 10^. The plot is similar to Fig. [ij the streamfunction extends into 
the domain where rapid distortion is important, and we therefore expect to see the shape 
of the streamfunction adjust to take account of the wave-turbulent interactions. As before, 
our computation of Zt enables us to develop an interpolation function X (z) in the gas: we 
take I (z) = 1 — e~(^/^^) , and for simplicity, we replace the a-dependent variable Zt (a) its 
average value, obtainable from Fig. fTOl 



B. Linear-stability analysis 

We carry out a stability analysis around the base state just constituted at a Reynolds 
number Re = 10^, and at an inverse Froude number Fr = 500. The inverse Weber number 
S is set to zero, a realistic assumption since the Froude number is large and thus gravity 
dominates over capillarity. The density and viscosity ratios are chosen such that the stability 
analysis models an air-water system under standard conditions: r = 1000 and m = 55 
respectively. These Reynolds and Froude numbers are chosen such that the so-called critical- 
layer instability is observed. We examine the growth rate of the disturbance with and 
without the PTS. A description of the change in the growth rate upon adding the PTS is 



a 


Amplitude of drift 


Not determined, although a value a — U^i/{2Uo) 
is given in the literature j25|. 


b 


Decay scale of the 
drift velocity 


Determined from a and the viscosity contrast 
m through the continuity of tangential stress 


C 


The kinetic energy amplitude 
is given by in wall units. 


C = 0.55, to agree with log-layer 
conditions in boundary-layer flow. 


n 


Exponent in the Van 
Driest damping function 


Chosen such that the Reynolds stress and the 
kinetic energy mimic wall turbulence as z ^ 0. 


A 


Length scale in the Van 
Driest damping function 


Chosen such that the linear region of the base 
flow is approximately 5 wall units in depth. 



Table II: Summary of the parameters used in Eqs. (32). 
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shown in Fig. [TTJ We also plot the wave speed in Fig. [TT} which is unchanged by the PTS- 
modelling. The segment of the dispersion curve that is of interest is the short-wave limit, 
for which a ^ 27r (where all lengths are measured relative to the gas-layer height do)- For 
longer waves, there is an interaction between the upper plate and the wave streamfunction, 
and the model no longer describes boundary-layer phenomena. The maximum growth rate 
occurs substantially above this lower cutoff value, at a wavenumber amax ^ 40, while gravity 
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Figure 11: The growth rate and wave speed of the PTS (continuous) and non-PTS waves (dotted 
lines), as a function of wavenumber a, normalized on the gas-layer thickness. 



stabilizes the wave above an upper critical wavenumber ac ^ 70. The location of the 
maximum growth rate is not changed by the PTS, although its value is shifted upwards, by 
about 10%. We do not continue the PTS dispersion curve below a = 10, since the PTS 
streamfunction is difficult to resolve numerically below this value. We do, however, describe 



the dispersion curve for the non-PTS case below this threshold in Fig. [12} There growth rate 
oscillates between positive and negative values. It is tempting to dismiss this as an effect of 
the upper plate; this is not the case however, since the effect persists upon increasing the 
depth of the gas layer. This anomalous region is far from the most dangerous mode, and we 
do not study it further. Now the structure of the flow field also changes as a consequence of 
the varying level of intensity of the rapid distortion, and it is to this variation that we now 
turn examining the perturbation velocity at different wavenumbers. 

a = 15: At this wavenumber, the wave parameters (growth rate and propagation speed) 

are 



Ar = 0.0004, 
Ar = 0.0008, 
Zr = 0.0065, 



0.2011, 
0.2011, 



with the PTS; 
with no PTS; 



(33) 



where Zc is the critical height, for which U (zc) = Cj.. The velocity field associated with this 
wave is shown in Fig. [T3| These figures all possess similar features, regardless of the PTS. 



The normal component of the velocity possesses successive extrema (maxima or minima). 
The first set of extrema is located at the surface z = 0, and, as in Sec. [IV| can be understood 
by examination of the kinematic condition, which in the moving frame of reference is simply 
dh/dt = w dit the interface. Upstream of a wave for which Cr > 0, the interfacial height 
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Figure 12: The growth rate at long wavelengths in the non-PTS case. 



decreases in time, implying that w < there; similarly, w > in the downstream region. In 
contrast to Sec. [TV] there is another set of extrema located above the critical layer; this is 
discussed below in the context of the a = 40 wave. In addition, since w = 
and u = ^ ^^ia{x-ct) ^ there is a 7r/2 phase difference between u and and thus the 

structure of the streamwise velocity can be understood simply as a phase shift relative to 
the normal velocity. The phase relationships discussed combine to give a distinctive phase 
relationship for the pressure field, shown in Fig. [M] In particular, at maximum growth, 
the pressure and the interfacial wave are approximately 7r/4 out of phase. Furthermore, the 
far-field structure of the PTS and non-PTS waves differ by the presence in the former case 
of successive velocity extrema, succeeding the one associated with the critical layer. These 
are due to the effects of the rapid distortion on the wave-induced motion, and are visible 



as velocity 'streaks' in the streamwise direction (see Fig. 13 (a)); these streaks are tilted 
forwards in the direction of the base fiow. 

a = 4:0: At this wavenumber, the wave parameters are 



A. 



0.0043, 



A, = 0.0047, 



= 0.0023; 



Cr = 0.1033, 
Cr = 0.1033, 



with the PTS; 
with no PTS; 



(34) 



We present an energy budget for the case a = 40 in Tab. Ill, which is calculated from the 



streamfunction using the balance law formulated in Eq. (26). The energy terms are further 



normalized such that the total energy of the PTS wave sums to unity. The two energy 



a = 40 


KINl 


KINg 


REYl 


REYg 


DISSl 


DISSg 


TURB 


NOR 


TAN 


PTS 


0.99 


0.01 


-0.04 


4.15 


-0.88 


-2.53 


0.04 


-2.06 


2.32 


No PTS 


0.99 


0.01 


-0.01 


2.81 


-0.62 


-1.76 





-1.01 


1.61 



Table III: The energy budget for the deep-water waves at o: = 40, near maximal growth. 



budgets possess only small relative differences, in spite of the large difference between the 
associated growth rates. This is because the kinetic energy possesses contributions both 
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Figure 13: Flow field structure for a — 15. The critical layer is shown as a dotted line above the 
interface. The effects of the rapid distortion are particularly visible in subfigures (a) and (b), which 
display the streamwise velocity and the streamwise velocity at x = 0, respectively. Subfigures (c) 
and (d) show the normal velocity and the normal velocity at x = 0. The effects of the rapid 
distortion are again visible in the pre-averaged version of the wave Reynolds stress, namely the 
product uw^ shown in (e) and (f). In each case, we have normalized the velocities by max \u\. 



from Cr and q, and while the difference in q is large, the difference in is small, and 
Cr 3> Q. Thus, the precise mechanism of excess instability is not visible from the energy 
budget. Nevertheless, the energy budget does indicate that the maximum destabilizing term 
is REYgj and we therefore study the wave Reynolds stress: 

n27T/a 
^0 

The sum REYq is thus 
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Figure 14: The pressure distribution for a — 15, (a), and a — 40, (b), normalized by max \u\. The 
critical layer is shown as a dotted line above the interface. 



and thus the wave Reynolds stress represents the energy associated with the power REYq. 



We plot the wave Reynolds stress in Fig. 15, The peak in r^ave next to the critical height Zc 




Figure 15: The unit-normalized wave Reynolds stress for a — 40. The critical layer is at Zc = 
0.0022, shown in the figure. The stress function for the non-PTS wave is given by the dotted 
line, while the continuous curve indicates the PTS wave. The PTS stress function possesses an 
oscillatory structure in the gas layer. 



represents a net transfer of energy from the mean flow U (z) into the perturbation flow. This 
is present in both the PTS wave and the non-PTS wave and indicates that the instability is a 
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critical-layer or Miles instability. Note, however that the PTS wave possesses an osciUation 
that propagates into the bulk gas flow, due to the interaction between the turbulence and 



the interface. In contrast to the a = 10 case, the velocity flelds in Fig. 16 possess the same 
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Figure 16: Flow field structure for a = 40. The critical layer is shown as a dotted line above the 
interface. The effects of the rapid distortion on the structure of the flow fleld no longer visible. The 
streamfunction decays in the bulk of the gas more rapidly than in the a = 10 case (as evidenced 



by a comparison between Figs. 16 (d) and [13| (d)). Subflgures (a) and (b) show the streamwise 
velocity and a proflle at x = 0; subflgures (c) and (d) show the normal velocity and a proflle at 
X = 0; subflgures (e) and (f) show the pre-averaged version of the wave Reynolds stress, namely 
the product uw. In each case, we have normalized the velocities by max \u\. 



structure, regardless of the PTS, and neither the oscillation nor the streamwise streaks are 



visible. This is consistent with the predictions of the piecewise-constant model in Sec. |IVC 
Finally, in Figs 



T7| and [18] we plot the turbulent kinetic energy and Reynolds stresses 
at a = 15 and 40 respectively. These are an order of magnitude smaller than the wave 
Reynolds stress uw^ conflrming the importance of the latter term, relative to the turbulent 
variables. Note that for the smaller of the two a- values, the maximum values of the turbulent 



variables lie close to the critical layer (Fig. 17), while for the larger a- value, the maximum 



lies far above the critical height (Fig. 18). This is consistent with the growth-rate curve in 
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Fig. 11, where the relative change in the growth rate (based on a comparison between the 



PTS no-PTS curves) was larger at smaller a-values, suggesting an interaction between the 
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Figure 17: Turbulence variables at o: = 15. The critical layer is shown as a dotted line above the 
interface, (a) and (b) show the turbulent kinetic energy; (c) and (d) show tangential Reynolds 
stress, while (e) and (f) show the normal Reynolds stress. We have normalized the perturbed 
turbulent quantities by max \u\, the maximum of the perturbed velocity. 



critical-layer mechanism and the wave turbulence. For larger a-values (in particular, close 
to the maximum growth rate), the relative change in the growth rate is smaller. 

In summary, the behaviour of the PTS wave differs from that of the non-PTS wave. The 
difference is both quantitative (the growth rate shifts), and qualitative (the flow field changes, 
especially at longer wavelengths). At longer wavelengths, the streamwise velocity exhibits 
streaks due to the interaction of the wave and the turbulence. At shorter wavelengths, this 
effect is reduced, although the growth rate is enhanced relative to the case without the PTS. 
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Figure 18: Turbulence variables si a — 40. The critical layer is shown as a dotted line above the 
interface, (a) and (b) show the turbulent kinetic energy; (c) and (d) show tangential Reynolds 
stress, while (e) and (f) show the normal Reynolds stress. We have normalized the perturbed 
turbulent quantities by max \u\^ the maximum of the perturbed velocity. 



C. Comparison with other work 

Since we have identified the instability as a critical-layer instability, it is appropriate to 
compare our results with the theory of Miles. Morland and Saffman [36j have developed 
explicit formulas that enable such a comparison. We also develop a comparison with deep- 
water waves based on the work of Boomkamp and Miesen in [18j. 

In the paper of Boomkamp and Miesen p!^J, the authors present one calculation of the 
wave Reynolds stress based on a boundary-layer gas profile and an exponential profile in the 
liquid. We perform this calculation again and construct a dispersion curve over a range of De- 
values. This result is shown in Fig. [191 and the wave Reynolds stress calculation is given as an 



34 







1 





OM 



1 




Figure 19: Comparison with the work of Boomkamp and Miesen [18j. The wave Reynolds stress 
at q:=40 is shown as in inset, with the critical layer at Zc = 0.0044. The growth-rate is designated 
by the solid curve, with a scale on the left-hand side, while the wave-speed is designated by the 
dashed-line curve, whose scale is on the right-hand side. 



inset. This calculation agrees with that given by Boomkamp and Miesen. More interesting is 
the dispersion curve. This is qualitatively similar to those given in Figs. [TT] and [T2| although 
the gravity, surface-tension and Reynolds numbers are different. Note in particular that the 
growth rate switches rapidly from positive to negative values at small a- values. This is not a 
finite-size effect, since it persists upon increasing the size of the computational domain and 
upon grid refinement. Having confirmed a qualitative similarity between the present work 
and a reconstruction of that of Boomkamp and Miesen, we endeavour to produce a more 
exact comparison. This we do by changing the parameter- values in the reconstructed work: 
specifically, we take Fr = 500, S = 0^ Re = 10^, and normalize the base-state velocity such 
that Ug {^) = |. This facilitates an accurate comparison between our work and that of 



Boomkamp and Miesen. We present the comparison in Fig. [20| together with a comparison 
with the Miles formula. 

Now the critical-layer theory of Miles [6] involves the solution of the Rayleigh equation 
with boundary conditions at infinity, and at a wavy, impermeable wall, which is supposed to 
represent the interface of two fluids with a large density contrast. For this model problem, 
Morland and Saffman [36] have derived an explicit but approximate formula for the growth 
rate in the case of the exponential profile 

where for comparison we take Uoo = and A = Re/Rel. Their approximate formula 

for the growth rate is 



rco (2 + aAf (4 + aAf \ f/oo 



= ^^^-TTj-r^-rri ( 1 - ) > co = J^. (35) 



We plot Ar against Cq for this simplified system, and compare the result with the relationship 
we have obtained between Ar and Cj.. The results are shown in Fig. [20| There is excellent 
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Figure 20: Comparison between our work, the analytical approximation of Morland and 
Saffman [36j, and a reconstruction of the work of Boomkamp and Miesen [18j. There is good 
agreement between all three models, although the analytical formula breaks down at small wave 
speeds, where surface tension exerts a stabilizing influence on the growth rate. 



agreement among the three curves shown, although the analytical curve corresponding to the 
formula (35) breaks down at small wave speeds, where surface tension stabilizes the system. 



Nevertheless, the form of curves in Fig. [20] is the same in each case, which, in addition to 
the energy budgets and Fig. [15} provides confirmation of the critical-layer nature of the 
instability. 

There are also some DNS results available in this field for comparison. In particular, we 
focus on the work of Sullivan et al and compare our results with DNS results found 
therein. While a direct comparison is not possible, since the work of Sullivan et al is for 
flow over a wavy wall, there our some qualitative similarities between the two two studies. 

Sullivan et al identify slow- wave and fast- wave cases. Thus, we make a comparison 
between our Fig. [7] and Fig. 17 in the work of Sullivan et a/., which is for slow waves, and for 
a critical layer that that plays no role in the dynamics. There is good qualitative agreement 
between these to sets of figures, although evidence of rapid distortion is absent in both 



cases. Next, we compare a fast-wave scenario, namely our Figs. [T3| and [T6| with Fig. 19 
in Sullivan et al. Again, there is good qualitative agreement between these sets of figures, 
in particular for the phase relationships at the interface (the corresponding results for the 
phase relationships of the pressure field also show good agreement). There is some evidence 
of rapid distortion in the secondary extrema in the streamwise velocity in Fig. 19 (a) of 
Sullivan et a/., although this is not conclusive: this comparison implies that our model over- 
estimates this effect. That said, these secondary oscillations cannot be reproduced by our 
linear theory if we neglect the PTS. 



D. Transition to the viscosity-contrast instability 

We have investigated two types of interfacial instability, and the effects of the PTS 
thereon. On the energy-budget side, these waves are distinguished by the energy term 
that produces the instability: either the viscosity-jump mechanism, or the critical-layer 
mechanism. On the turbulence side, they are distinguished by the wave speed: the viscosity- 
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jump waves are slow, while the critical-layer waves are fast. This separation of speeds leads 
to distinct properties with respect to the PTS. We want to find a transition from one 
regime to another: this can be achieved by changing the gravity number Fr (inverse Froude 
number). We study the deep-water waves again. We obtain the dispersion curve for a range 
of gravity numbers, and at fixed Reynolds number, and reduce the gravity number. This can 
be achieved either by changing the degree of density stratification, the gas-layer thickness, 
or by changing the shear velocity Uq: 



Fr = 



9 {pL - Pg) dp 
PgUI ' 



The dispersion curves are shown as a function of gravity number, expressed as a fraction of 
Fro = 500 in Fig. 21, The growth rate changes shape as the gravity number is decreased: 



the high-gravity number shape is similar to that observed in this section for critical-layer 



waves, while the low- gravity number shape is similar to that observed in Sec. IV for viscosity 
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Figure 21: The growth rate and wave speed in the absence of the PTS, as a function of gravity 
number, for Ftq = 500 and Re = 10^. The instability changes character as the gravity number is 
decreased. 



stratified waves. The plot of wave speed in Fig. 21 (b) indicates that the high-gravity number 
waves are fast while the low-gravity number waves are slow. Going to larger Fr-values 
stabilizes the interface completely, although we have not shown this effect in Fig. 21 For 
each dispersion curve, we obtain the energy budget associated with maximum growth, given 
in Tab. [IV| We observe a transition from critical-layer to viscosity-driven instability, with 
decreasing gravity number. The term TAN represents only a small positive contribution 
at high Fr-values, while at low Fr-values, it is the dominant positive contribution. The 
presence of the PTS does little to modify this partition of energy between the tangential 
and REYg terms. 

This mechanism for modifying the character of the instability could be applied to the 
thin-film case, with one reservation. To move the critical layer sufficiently far into the bulk 
gas domain such that U" {z^) is significant, it is necessary to increase the Reynolds number. 
This in turn will destabilize the liquid through an internal mode, which necessitates the 
turbulent modelling both of the base liquid fiow, and of the liquid PTS, which is beyond 
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the scope of the present work. However, by naively retaining the hnear profile in the liquid, 
we have observed a switch between critical-layer and viscosity-stratified waves in the thin- 
film case, through a suitable modification of the gravity and Reynolds numbers {Re = 10^, 
Fr = 0.1, maximum growth rate at a = 25). 



Fr 


REYl 


REYg 


DISSl 


DISSg 


TURB 


NOR 


INT 


Fro 


0.00 


1.75 


-0.39 


-1.10 


0.00 


-0.58 


1.00 


\Fro 


-0.01 


0.42 


-0.32 


-0.74 


0.00 


-0.17 


1.00 


iFro 


0.00 


0.07 


-0.19 


-0.72 


0.00 


-0.04 


1.00 




0.00 
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0.00 
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0.00 
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NOR 


INT 


Fro 
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\Fro 
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0.00 
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Table IV: Energy budget detailing the transition from critical-layer to viscosity-stratified waves, as 
a function of gravity number, where Ftq = 500 and Re = 10^. The budgets have been normalized 
such that TAN=1 in each case. In the first table, we have included the PTS; in the second table, 
they are set to zero. 



VI. CONCLUSIONS 

We have investigated the stability of an interface separating a liquid layer from a fully- 
developed turbulent gas fiow. The linear-stability analysis involves the study of the dynamics 
of a wave on the interface, and this wave interacts with the turbulence and induces pertur- 
bation turbulent stresses (PTS), which modify the stability properties of the system. Using 
a separation-of-domains technique, we derived a model for the PTS, based on the Orr~ 
Sommerfeld equation for the streamfunction. We also developed a model of the base fiow 
that takes near- wall (interface) regions into account, and provides the friction velocity U^i 
as a function of Reynolds number. 

We have applied our model in two distinct cases. For fiow over a thin viscous film, and 
at moderate values of the Reynolds number, the turbulence modelling does not materially 
affect the growth rate, although the structure of the velocity field is modified: the streamwise 
velocity develops streaks that extend into the bulk gas layer, as observed in DNS On the 
other hand, for deep-water waves and at high Reynolds numbers, the maximum growth rate 
is shifted upwards by the PTS, and the fiow structure is again modified, especially at longer 
wavelengths, when the spatial extent of the streamfunction extends into the rapid-distortion 
domain. The waves in the thin film and the deep channel differ are slow and fast respectively, 
compared with the shear velocity at the upper plate. They can also be classified respectively 
as viscosity-stratified, or critical-layer waves. The waves observed in both cases can, however. 
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be brought into co-incidence by a modification of the gravity number. By decreasing the 
gravity number in the deep-water waves, we have observed a transition from the critical-layer 
to the viscosity-stratified waves. The transition has, effectively, been effected by modifying 
the shear velocity, the gas-layer depth, or the degree of stratification. This suggests that 
a detailed parameter study will be useful in understanding the different mechanisms that 
generate instability, an approach we develop in O Naraigh et al |[3j 
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